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Abstract. We investigate the strong stability preserving (SSP) property of two-step Runge— 
Kutta (TSRK) methods. We prove that all SSP TSRK methods belong to a particularly simple 
subclass of TSRK methods, in which stages from the previous step are not used. We derive simple 
order conditions for this subclass. Whereas explicit SSP Runge-Kutta methods have order at most 
four, we prove that explicit SSP TSRK methods have order at most eight. We present TSRK 
methods of up to eighth order that were found by numerical search. These methods have larger SSP 
coefficients than any known methods of the same order of accuracy, and may be implemented in a form 
with relatively modest storage requirements. The usefulness of the TSRK methods is demonstrated 
through numerical examples, including integration of very high order WENO discretizations. 

1. Strong Stability Preserving Methods. The concept of strong stability 
preserving methods was first introduced by Shu and Osher in [40] for use with total 
variation diminishing spatial discretizations of a hyperbolic conservation law: 

u t + f(U) x = 0. 

When the spatial derivative is discretized, we obtain the system of ODEs 

u t =F(u), (1.1) 

where u is a vector of approximations to U, uj w U(xj). The spatial discretization 
is carefully designed so that when this ODE is fully discretized using the forward 
Euler method, certain convex functional properties (such as the total variation) of 
the numerical solution do not increase, 

\\u n + AtF(u n )\\ < \\u n \\ (1.2) 

for all small enough step sizes At < AipE- Typically, we need methods of higher order 
and we wish to guarantee that the higher-order time discretizations will preserve this 
strong stability property. This guarantee is obtained by observing that if a time 
discretization can be decomposed into convex combinations of forward Euler steps, 
then any convex functional property (referred to herein as a strong stability property) 
satisfied by forward Euler will be preserved by the higher-order time discretizations, 
perhaps under a different time-step restriction. 

Given a semi-discretization of the form (1.1) and convex functional || • ||, we assume 
that there exists a value A£fe such that, for all u, 

\\u + AtF(u)\\ < \\u\\ for < At < At FE . (1.3) 
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A /c-step numerical method for (1.1) computes the next solution value u n+1 from 
previous values u n ~ k+1 , . . . , u n . We say that the method is strong stability preserving 
(SSP) if (in the solution of (1.1)) it holds that 

<m^{\\u n \\,\\u n - x \\,...,\\u n - k+l \\}. (i.4) 

whenever (1.3) holds and the timestep satisfies 

At < CAt FE . (1.5) 

Throughout this work, C is taken to be the largest value such that (1.5) and (1.3) 
together always imply (1.4). This value C is called the SSP coefficient of the method. 
For example, consider explicit multistep methods [39]: 

k 

u n+1 = (a i u n+1 - i + Atp,F{u n+1 -' 1 )) . (1.6) 

i=l 

l 

Since 2_/j=i a i = 1 f° r an y consistent method, any such method can be written as 
convex combinations of forward Euler steps if all the coefficients are non-negative: 

k 



u 

»=l 



"+ 1 = Y on ( u n+1 ~* + — AiF(w n+1_i ) 
^ v at 



If the forward Euler method applied to (1.1) is strongly stable under the timestep 
restriction At < AiFE and ai, Pi > then the solution obtained by the multistep 
method (1.6) satisfies the strong stability bound (1.4) under the timestep restriction 

At < min — A^fe , 

i Pi 

(if any of the /3's are equal to zero, the corresponding ratios are considered infinite). 
In the case of a one-step method the monotonicity requirement (1.4) reduces to 

IK +1 ||<IKI|. 

For example, an s-stage explicit Runge-Kutta method is written in the form [40], 



1—1 

Y,(a ijU W + AtfajFiu^j) , (1.7) 



u n+1 =u^. 

If all the coefficients are non-negative, each stage of the Runge-Kutta method can 
be rearranged into convex combinations of forward Euler steps, with a modified step 
size: 

i-l 

\\u^\\ = \\J2(^^+At^F(u^)) || 

i-l 

3=0 



U U) + At^F(u^) 



Now, since each \\u^ + At^F(u^)\\ < as long as ^At < At FE , and since 

^2]=a a ij — 1 by consistency, we have < as long as ^-At < Ai FE for 

all i and j. Thus, if the forward Euler method applied to (1.1) is strongly stable 
under the timestep restriction At < AtpE, i.e. (1.3) holds, and if a,j, fiij > then 
the solution obtained by the Runge-Kutta method (1.7) satisfies the strong stability 
bound (1-2) under the timestep restriction 

At < min^-At FE . 

»>i Pij 

As above, if any of the /3's are equal to zero, the corresponding ratios are considered 
infinite. 

This approach can easily be generalized to implicit Runge-Kutta methods and 
implicit linear multistep methods. Thus it provides sufficient conditions for strong 
stability of high-order explicit and implicit Runge-Kutta and multistep methods. In 
fact, it can be shown from the connections between SSP theory and contractivity 
theory [11, 12, 19, 20] that these conditions are not only sufficient, they are necessary 
as well. 

Research in the field of SSP methods focuses on finding high-order time discretiza- 
tions with the largest allowable time-step. Unfortunately, explicit SSP Runge-Kutta 
methods with positive coefficients cannot be more than fourth-order accurate [32, 38], 
and explicit SSP linear multistep methods of high-order accuracy require very many 
steps in order to have reasonable timestep restrictions. For instance, to obtain a fifth- 
order explicit linear multistep method with a time-step restriction of At < 0.2AipE 
requires nine steps; for a sixth-order method, this increases to thirteen steps [34]. In 
practice, the large storage requirements of these methods make them unsuitable for 
the solution of the large systems of ODEs resulting from semi-discretization of a PDE. 
Multistep methods with larger SSP coefficients and fewer stages have been obtained 
by considering special starting procedures [22, 37]. 

Because of the lack of practical explicit SSP methods of very high order, high- 
order spatial discretizations for hyperbolic PDEs are often paired with lower-order 
time discretizations; some examples of this include [5, 6, 7, 9, 10, 27, 33, 36, 43]. 
This may lead to loss of accuracy, particularly for long time simulations. In an 
extreme case [13], WENO schemes of up to 17th-order were paired with third-order 
SSP Runge-Kutta time integration; of course, convergence tests indicated only third- 
order convergence for the fully discrete schemes. Practical higher-order accurate SSP 
time discretization methods are needed for the time evolution of ODEs resulting from 
high-order spatial discretizations. 

To obtain higher-order explicit SSP time discretizations, methods that include 
both multiple steps and multiple stages have been considered. These methods are 
a subclass of explicit general linear methods that allow higher order with positive 
SSP coefficients. Gottlieb et. al. considered a class of two-step, two-stage methods 
[14]. Another class of such methods was considered by Spijker [41]. Huang [21] 
considered hybrid methods with many steps, and found methods of up to seventh- 
order (with seven steps) with reasonable SSP coefficients. Constantinescu and Sandu 
[8] considered two- and three-step Runge-Kutta methods, with a focus on finding SSP 
methods with stage order up to four. 

In this work we consider a class of two-step multi-stage Runge-Kutta methods, 
which are a generalization of both linear multistep methods and Runge-Kutta meth- 
ods. We have found that deriving the order conditions using a generalization of 
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the approach presented in [2], and formulating the optimization problem using the 
approach from [31] allows us to efficiently find methods of up to eighth order with 
relatively modest storage requirements and large effective SSP coefficient. We also 
report optimal lower-order methods; our results agree with those of [8] for second, 
third, and fourth-order methods of up to four stages, and improve upon other meth- 
ods previously found both in terms of order and the size of the SSP coefficient. 

The major result of this paper is the development of SSP two-step Runge-Kutta 
methods of up to eighth order that are efficient and practical. In Section 2, we 
discuss some classes of two-step Runge-Kutta (TSRK) methods and prove that all 
SSP TSRK methods belong to one of two simple subclasses. In Section 3, we derive 
order conditions and show that explicit SSP TSRK methods have order at most eight. 
In Section 4, we formulate the optimization problem, give an efficient form for the 
implementation of SSP two-step Runge-Kutta methods, and present optimal methods 
of up to eighth order. The properties of our methods are compared with those of 
existing SSP methods including Runge-Kutta, linear multi-step, and hybrid methods 
[21], as well as the two- and three-step methods in [8]. Numerical verification of the 
optimal methods and a demonstration of the need for high-order time discretizations 
for use with high-order spatial discretizations is presented in Section 5. Conclusions 
and future work are discussed in Section 6. 

2. SSP Two-step Runge Kutta Methods. The principal focus of this work 
is on the strong stability preserving properties of two-step Runge-Kutta (TSRK) 
methods. A general class of TSRK methods was studied in [25, 4, 16, 45]. TSRK 
methods are a generalization of Runge-Kutta methods that include values and stages 
from the previous step: 

s s 

y? = diU 71 - 1 + (1 - di)u n + At ^FWf 1 ) + Ai^ayify?), 1 < * < *, 

3=1 3=1 

(2.1a) 

S S 

u n+i _ 9u n-i + 0j u „ + Ai^Sji^- 1 ) + At^bjFiy^). (2.1b) 

3=1 3=1 



Here u n and u n 1 denote solution values at the times t = nAt and t = (n — l)At, 
while the values yf are intermediate stages used to compute the solution at the next 
time step. We will use the matrices and vectors A, A, b, b, and d to refer to the 
coefficients of the method. 

We are interested only in TSRK methods that have the strong stability preserving 
property. As we will prove in Theorem 1, this greatly reduces the set of methods 
relevant to our study. Except in special cases, the method (2.1) cannot be strong 
stability preserving unless all of the coefficients a>ij,bj are identically zero. A brief 
explanation of this requirement is as follows. Since method (2.1) does not include 
terms of the form y™ -1 , it is not possible to write a stage of method (2.1) as a convex 
combination of forward Euler steps if the stage includes terms of the form F(y"^ 1 ). 
This is because those stages depend on u n ~ 2 , which is not available in a two-step 
method. 

Hence we are led to consider simpler methods of the following form (compare [15, 
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p. 362]). We call these Type I methods: 

s 

d^- 1 + (l-di)u n + At^a %0 F{yJ), l<i<s, (2.2a) 

3=1 

S 

On 71 - 1 + (1 - 8)u n + AtJ2 b J F (yj)- ( 2 - 2b ) 

3=1 

Now consider the special case in which the method (2.1) has some stage yf iden- 
tically equal to u n . Then including terms proportional to F{u n ) will not prevent the 
method from being written as a convex combination of forward Euler steps; further- 
more, since y™" 1 = u™ -1 , terms of the form F(u n ~ 1 ) can also be included. This leads 
to what we will call Type II methods, which have the form: 

u n , (2.3a) 

S 

din 71 ' 1 + (1 - di)u n + a i AtF(u n - 1 ) + AtJ2 a >ij F (yj)> 2 < * < s > ( 2 - 3b ) 

3=1 

S 

du n - 1 + (1 - 6)u n + biAtF(u n ^ 1 ) + AtJ2 b i F (Vj)- (2.3c) 

3=1 

Here we have assumed that the first stage is the one equal to u n , which involves no loss 
of generality. We can refer to the coefficients of Type II methods in the matrix/vector 
notation of (2.1) except that matrix A reduces to vector a and we have d% = a\ = 
and aij = for all 1 < j < s. 

Remark 1. From a theoretical perspective, the distinction between Type I and 
Type II methods may seem artificial, since the class of all Type I methods is equiv- 
alent to the class of all Type II methods. From a practical perspective, however, the 
distinction is very useful. Transforming a given method from one type to the other 
generally requires adding a stage. Thus the class of s-stage Type I methods and the 
class of s-stage Type II methods are distinct (though not disjoint). So it is natural to 
refer to a method as being of Type I or Type II, depending on which representation 
uses fewer stages; this convention is used throughout the present work. 

Remark 2. Type I methods (2.2) and Type II methods (2.3) are equivalent to 
the (two-step) methods of Type 4 and Type 5, respectively, considered in [8]. 

2.1. The Spijker Form for General Linear Methods. TSRK methods are 
a subclass of general linear methods. In this section, we review the theory of strong 
stability preservation for general linear methods [41]. A general linear method can be 
written in the form 

/ m 

• A/ )J /,,/•: „■;'). (l<*<m), (2.4a) 

3 = 1 3=1 

w n j 3 , (i < i < 0- ( 2 - 4b ) 

The terms x™ are the I input values available from previous steps, while the wj includes 
both the output values and intermediate stages used to compute them. Equation 
(2.4b) indicates which of these values are used as inputs in the next step. 

We will frequently write the coefficients Sij and Uj as a m x I matrix S and a 
m x m matrix T, respectively. Without loss of generality (see [41, Section 2.1.1]) we 
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y? = 

u n+1 = 



Vi = 

v n i = 

u n+1 = 



-n+1 



assume that 



Se = e. 



(2.5) 



where e is a vector with all entries equal to unity. This implies that every stage is a 
consistent approximation to the solution at some time. 

Runge-Kutta methods, multi-step methods, and multi-step Runge-Kutta meth- 
ods are all subclasses of general linear methods, and can be written in the form (2.4). 
For example, an s-stage Runge-Kutta method with Butcher coefficients A and b can 
be written in form (2.4) by taking I = 1, m = s + 1, J = {m}, and 

S = (1,1,...,1) T , T= ( * ° 

Linear multistep methods 



3=1 



3=0 



admit the Spijker form 
(1 


... 

V"; a i-i 



'•■ 
1 
. .. aij 



T (;+i)x(;+i) _ 



/0 

o '■• 

... 



■. 



... /v 



where I is the number of steps, m = I + 1, and J = {2, ...,/ + 1}. 

General TSRK methods (2.1) can be written in Spijker form as follows: set m 
2s + 2, / = s + 2, J = {s + 1, s + 2, . . . , 2s + 2}, and 

>T 



— l n ,n—l 



x" = (u- ',y{ 

(0 1 \ 

1 

d e d 



(2.6a) 
(2.6b) 



/ 0\ 




\e o 1-6) 



A A 

\b T o b T o/ 



(2.6c) 



Type I methods (2.2) can be written in a simpler form with m — s + 2, I 
J = {l,s + 2}, and 





(«- 




w" = (u n 


Vi 


2/2™, ■ 









(° 







d 


e 




T = 


A 


:) 





1 




\o 


b T 





Type II methods (2.3) can also be written in a simpler form with m = s 
I = 2, J= {2,s + 2}: 




lU "-V n ,2/£,...,yI> n+1 ) 
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2.2. The SSP Coefficient for General Linear Methods. In order to analyze 
the SSP property of a general linear method (2.4), we first define the vector f = 
[F(wi), F(w2), ■ • ■ , F(w m )] T , so that (2.4a) can be written compactly as 

w = Sx + AiTf. (2.7) 

Adding rTw to both sides of (2.7) gives 



(I + rT) w = Sx + rT ( w + — f 



r 



Assuming that the matrix on the left is invertible we obtain, 

w = (I + rT) _1 Sx + r(I + rT) _1 T + — f 

= Rx + P (w + — f J , (2.8) 



r 



where we have defined 



P = r(I + rT)- lr r, R= (I + rT^S = (I-P)S. (2.9) 

Observe that, by the consistency condition (2.5), the row sums of [R P] are each 
equal to one: 

Re + Pe = (I P)Se + Pe = e Pe + Pe = e. 

Thus, if R and P have no negative entries, each stage Wi is given by a convex com- 
bination of the inputs Xj and the quantities Wj + (At/r)F(wj). In other words, this 
method is a convex combination of forward Euler steps. Hence any strong stability 
property of the forward Euler method is preserved by the method (2.7) under the 
time step restriction given by At < C(S, T)AtpE where C(S, T) is defined as 

C(S, T) = sup {r:(I + rT)" 1 exists and P > 0, R > 0} , 

r 

where P and R are defined in (2.9). By the foregoing observation, it is clear that the 
SSP coefficient of method (2.8) is greater than or equal to C(S, T). 

To state precisely the conditions under which the SSP coefficient is, in fact, equal 
to C(S, T), we must introduce the concept of reducibility. A Runge-Kutta method is 
said to be reducible if there exists a method with fewer stages that always produces 
the same output. One kind of reducibility is known as US -reducibility; a Runge-Kutta 
method is HS-reducible if two of its stages are identically equal. This definition is 
extended in a natural way to general linear methods in [41, Theorem 3.1]; hereafter 
we refer to the reducibility concept defined there also as HS -reducibility. 

Lemma 1. ([41, Theorem 3.1]) Let S,T be an US -irreducible representation of a 
general linear method. Then the SSP coefficient of the method is C = C(S, T). 

2.3. Restrictions on the coefficients of SSP TSRK methods. In light of 
Lemma 1, we are interested in methods with C(S,T) > 0. The following lemma 
characterizes such methods. 
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Lemma 2. ([41, Theorem 2.2(i)]) C(S,T) > if and only if all of the following 
hold: 



S > 0, 

T > 0, 
Inc(TS) < Inc(S), 
Inc(T 2 ) < Inc(T). 



(2.10a) 
(2.10b) 
(2.10c) 
(2.10d) 



where all the inequalities are element-wise and the incidence matrix of a matrix M 
with entries rriij is 

Inc(M)ij- = < J l]T 

I if rriij = 0. 



To apply Lemma 1, it is necessary to write a TSRK method in HS-irreducible 
form. A trivial type of HS-reducibility is the case where two stages y™, y™ are identi- 
cally equal; i.e., where the following condition holds for some i ^ j: 



di = dj, rows i,j of A arc identical, and rows i,j of A are identical (2-11) 

This type of reducibility can be dealt with by simply combining the two stages; hence 
in the following theorem we assume any such stages have been eliminated already. 
Combining Lemma 1 and Lemma 2, we find that all SSP TSRK methods can be 
represented as Type I and Type II methods, introduced in Section 1. 

Theorem 1. Let S,T be the coefficients of a s-stage TSRK method (2.1) in the 
form (2.6) with positive SSP coefficient C > such that (2.11) does not hold for any 
i ^ j. Then the method can be written as an s-stage HS-irreducible method of Type I 
or Type II. 

Proof. Consider a method satisfying the stated assumptions. Examination of S, T 
reveals that the method is HS-irreducible in form (2.6) if there is no yj identically 
equal to u n . In this case, we can apply Lemma 1 to obtain that C(S,T) > 0. Then 
condition (2.10c) of Lemma 2 implies that A = b = 0. Under this restriction, methods 
of the form (2.1) simplify to Type I methods (2.2). 

Now consider the case in which yj = u n for some j. If necessary, reorder the 
stages so that y™ = u n . Then rows s + 1 and s + 2 of [S T] in the representation (2.6) 
are equal. Thus we can rewrite the method in form (noting that also y™ _1 = u™ -1 ) 
as follows: Set m = 2s + 1, I = s + 1, J = {s, s + 1, s + 2, . . . , 2s + 2}, and 



x" 
w" = ( 



■ii. 



7i—l ,,n\ 



n—1 n n — l .,71 — 1 



fl \ 

01 

d e d 

\e 1-0/ 



n — 1 ni n „,n 



,2/2 ,2/3 ,--->2/ s ,2/1,2/2, ■••,!/,! 







Ai A 2:s 










A 

b T 0/ 



(2.12a) 
(2.12b) 

(2.12c) 



Here Ai and b\ represent the first column and first element of A and b, respectively, 
while A2 :s and bj. s represent the remaining columns and remaining entries. Since 
(2.12) is HS-irreducible, we can apply Lemma 1. Then we have that C(S,T) > 0, so 



that Lemma 2 applies. Applying condition (2.10c) of Lemma 2 to the representation 
(2.12), we find that A.2- s and bj s must vanish, but Ai and b\ may be non-zero. The 
resulting methods are HS-irreducible Type II methods (2.3). □ 

Remark 3. One could formulate a version of Theorem 1 with the hypothesis that 
the method in form (2.1) is HS-irreducible. This would lead only to the class of Type 
I methods, and in the resulting formalism the number of stages for Type II methods 
would be artificially increased by one (see Remark 1). The advantage of explicitly 
considering the case y™ — u n is that we obtain Type II methods with the number of 
stages that accurately represents their cost. 

We now introduce a compact, unified notation for Type I and Type II methods. 
This form is convenient for expressing the order conditions and restrictions on the 
coefficients. First we rewrite an s-stage Type II method (2.3) by including u™" 1 as 
one of the stages: 

v n i = u n , 

s 

y? = d^- 1 + (1 - di)u n + Ai^OtfFfo?), 2 < i < s, 



ir ' = BvJ 1 - 1 + (1 - 6)u n + AtJ^&jFfo?). 



3=0 

Then both Type I and Type II methods can be written in the compact form 

y n = du™- 1 + (1 - d)u n + AtAr, (2.13a) 

u n + l _ &u n-l Qy, + Ai b T f«, (2.13b) 

where, for Type I methods the coefficients with bars are equal to the corresponding 
coefficients without bars in (2.2) and 

y" = [y?,...,y s "] T , f n = [F(v?),... 1 FW)] T . 

Meanwhile, for Type II methods 

y" = [u n -\u n , y 2 ", . . . , V :\ T , f" = [F^"- 1 ), F(u n ),F(y^), . . . , F(y^)] T , 
d= [l,0,d 2 ,...,4] T , b=[6xb T ] T , A-'° ° 



a a; ' 



where dj, b, A, b\, a refer to the coefficients in (2.3). 

It is known that irreducible strong stability preserving Runge-Kutta methods 
have positive stage coefficients, > and strictly positive weights, bj > 0. The 
following theorem shows that similar properties hold for SSP TSRK methods. The 
theorem and its proof are very similar to [32, Theorem 4.2]. In the proof, we will use 
a second irreducibility concept. A method is said to be DJ-reducible if it involves one 
or more stages whose value does not affect the output. If a TSRK method is neither 
HS-reducible nor DJ-reducible, we say it is irreducible. 

Theorem 2. The coefficients of an HS-irreducible TSRK method of Type I (2.2) 
or Type II (2.3) with positive SSP coefficient satisfy the following bounds: 

A > 0, b > 0, < d < 1, and < 6 < 1. (2.14) 
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Furthermore, if the method is also D J -irreducible, the weights must be strictly positive: 



b > 0. 



(2.15) 



All of these inequalities should be interpreted component-wise. 

Proof. Application of Lemma 1 implies that C(S, T) > 0. Therefore Lemma 2 
applies. The first result (2.14) then follows from conditions (2.10a) and (2.10b) of 
Lemma 2. To prove the second part, observe that condition (2.10d) of Lemma 2 
means that if bj — for some j then 



Since A,b are non-negative, (2.16) implies that either bi or ay is zero for each value 
of i. Now partition the set S = {1, 2, . . . , s} into Si, S2 such that bj > for all j 6 Si 
and bj = for all j G 1S2. Then ay = for all i G Si and j G £2- This implies that 
the method is DJ-reducible, unless 6>2 = 0. 

3. Order Conditions and a Barrier. Order conditions for TSRK methods 
up to order 6 have previously been derived in [25]. However, two of the sixth-order 
conditions therein appear to contain errors (they do not make sense dimensionally). 
Alternative approaches to order conditions for TSRK methods, using trees and B- 
series, have also been identified [4, 16]. 

In this section we derive order conditions for TSRK methods of Types I and II. 
The order conditions derived here are not valid for for the general class of methods 
given by (2.1). Our derivation follows Albrecht's approach [2], and leads to very 
simple conditions, which are almost identical in appearance to order conditions for 
RK methods. For simplicity of notation, we consider a scalar ODE only. For more 
details and justification of this approach for systems, see [2]. 

3.1. Derivation of Order Conditions. When applied to the trivial ODE 
u'(t) = 0, any TSRK scheme reduces to the recurrence u n+1 = du n ~ x + (1 — 6)u n . 
For a SSP TSRK scheme, we have < 6 < 1 (by Theorem 2) and it follows that the 
method is zero-stable. Hence to prove convergence of order p, it is sufficient to prove 
consistency of order p (see, e.g., [24, Theorem 2.3.4]). 

Let u(t) denote the exact solution at time t and define 



where c identifies the abscissae of the TSRK scheme. These represent the correct 
stage values and the corresponding correct values of F. Then the truncation error r™ 
and stage truncation errors r" are implicitly defined by 




(2.16) 



y" = [u(t n + ciAt), . . . , u(t n + c s At)] 
f" = [F(u(t n + Cl Atj), F{u{t n + c s At))}, 



u 



y" = du"- 1 + (e - d)u n + AtAi n + Atr n , 
(t n+ i) = On 11 - 1 + (1 - 8)u n + Atb T i n + Air™. 



(3.1a) 
(3.1b) 
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To find formulas for the truncation errors, we make use of the Taylor expansions 

oo 

u(t n + Cl At) = J2 77At fe C^W(t„), 



k\ 



oo 

F(u(t n + aAt)) = u'(t n + c t At) = (fc _ 1) , At fc - 1 c^ 1 ^ fc )(t n ), 



k=Q 



J(fc-l)! 
At 

fcT 



Substitution gives 



r n = ^r fe At fc - 1 uW(t n ), (3.2a) 
fe=i 

oo 

r n =Y,r k At k - 1 u^(t n ), (3.2b) 



fc=i 



where 



-St 1 -!-"'"-^" 

Subtracting (3.1) from (2.13) gives 

e" = de"- 1 + (e - d)e" + AtAS n - Air", (3.3a) 
e n+i _ 9f n-i + ^ _ fl j e „ + At ^s n - Atr n , (3.3b) 

where e n+1 = u n+1 — u(t n +i) is the global error, e n = y™ — y™, is the global stage 
error, and S n = f™ — f" is the right-hand-side stage error. 

If we assume an expansion for the right-hand-side stage errors S n as a power series 
in At 

p-i 

S n = S k Atk + 0(At p ), (3.4) 

k=Q 

then substituting the expansions (3.4) and (3.2) into the global error formula (3.3) 
yields 

p-i P 
e" = de"- 1 + (e - d)e" + ^ A^Ai fe+1 - ^ T k u {k) (t n )At k + 0(At p+1 ), 

k=0 k=l 

(3.5a) 

p-i p 

e n + l = 0e n-l + (1 _ Q yn + b T d%At k+1 - ^ T k U^ (t n ) At k + 0{At p+1 ). 

k=0 k=l 

(3.5b) 

Hence we find the method is consistent of order p if 

h T S k =0 (0<k<p-l) and r k = (1 < k < p). (3.6) 
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It remains to determine the vectors in the expansion (3.4). In fact, we can 
relate these recursively to the e&. First we define 

t n = t„e + cAt, 
F(y,t) = [F(y 1 {t 1 )),...,F(y s (t s ))] T . 

Then we have the Taylor series 

oo 1 

r = F( y ",t„) = r + £ -(y n - y n y • F (j) (y n ,t„) 

3=1 J ' 
oo 1 

3=1 ■'' 

where 

F^(y,t) = [F^(y 1 (t 1 )),...,F^(y s (t s ))] T , 
Ej (t) = [F^(y(t 1 )),...,F^(y(t s ))} T , 
and the dot product denotes component-wise multiplication. Thus 

5 n = f „ _ f „ = J- -(e n y ■ Sj (t n e + cAt). 

Since 

gj (^e + c)=£— C<gf (i„), 
where C = diag(c), we finally obtain the desired expansion: 

* n = EEw c,(en)i,g ? )(<n) - (3J) 

To determine the coefficients <5fc, we alternate recursively between (3.7) and (3.5a). 
Typically, the abscissae c are chosen as c = Ae so that T\ = 0. With these choices, 
we collect the terms relevant for up to fifth-order accuracy: 



Terms 


appearing 


in #i 





Terms 


appearing 


in e 2 




Terms 


appearing 


in 82 


T2 


Terms 


appearing 


in e 3 


At 2 , T 3 


Terms 


appearing 


in £ 3 


Ct 2 ,At 2 ,t 3 _ 


Terms 


appearing 


in 64 


ACr 2 , A 2 t 2 , At 3 ,t 4 


Terms 


appearing 


in £4 


ACt 2 , A 2 t 2 , At 3 , t 4 , CAr 2 , Cr 3 , C 2 r 2 , r 2 . 



The order conditions are then given by (3.6). In fact, we are left with order 
conditions identical to those for Runge-Kutta methods, except that the definitions of 
the stage truncation errors Tfe,Tfc, and of the abscissas c are modified. For a list of 
the order conditions up to eighth order, see [1, Appendix A]. 
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3.2. Order and Stage Order of TSRK Methods. The presence of the term 
r| in <$4 leads to the order condition b T r| = 0. For SSP methods, since b > (by 
Theorem 2), this implies that t\ = 0, i.e. fifth-order SSP TSRK methods must have 
stage order of at least two. The corresponding condition for Runge-Kutta methods 
leads to the well-known result that no explicit RK method can have order greater than 
four and a positive SSP coefficient. Similarly, the conditions for seventh order will 
include b T r§ = 0, which leads (together with the non-negativity of A) to the result 
that implicit SSP RK methods have order at most six. In general, the conditions 
for order 2k + 1 will include the condition b T r^ = 0. Thus, like SSP Runge-Kutta 
methods, SSP TSRK methods have a lower bound on the stage order, and an upper 
bound on the overall order. 

Theorem 3. Any TSRK method (2.1) of order p with positive SSP coefficient 
has stage order at least L^ipJ- 

Proof. Without loss of generality, we consider only irreducible methods; thus we 
can apply Theorems 1 and 2. Following the procedure outlined above, we find that 
for order p, the coefficients must satisfy 

h T T 2 k = Q, fc = l,2,..., 

Since b > by Theorem 2, this implies that 

T 2 k = 0, k=l,2,..., 

□ 

Application of Theorem 3 dramatically simplifies the order conditions for high 
order SSP TSRKs. This is because increased stage order leads to the vanishing of 
many of the order conditions. Additionally, Theorem 3 leads to an upper bound on 
the order of explicit SSP TSRKs. 

Theorem 4. The order of an explicit SSP TSRK method is at most eight. Fur- 
thermore, if the method has order greater than six, it is of Type II. 

Proof. Without loss of generality, we consider only irreducible methods; thus we 
can apply Theorems 2 and 3. 

To prove the second part, consider an explicit irreducible TSRK method with 
order greater than six. By Theorem 3, this method must have stage order at least 
three. Solving the conditions for stage yi to have stage order three gives that C2 must 
be equal to —1 or 0. Taking c-i — — 1 implies that ij2 = u n_1 , so the method is of type 
II. Taking C2 = implies y2 = y\ = u n ; in this case, there must be some stage yj not 
equal to u n and we find that necessarily Cj = — 1 and hence yj = u n ~ l . 

To prove the first part, suppose there exists an irreducible SSP TSRK method 
(2.13) of order nine. By Theorem 3, this method must have stage order at least four. 
Let j be the index of the first stage that is not identically equal to u™ -1 or u n . Solving 
the conditions for stage yj to have order four reveals that Cj must be equal to —1,0, 
or 1. The cases Cj — —1 and Cj — lead to yj = and yj = u n , contradicting 

our assumption. Taking Cj — 1 leads to dj ~ 5. By Theorem 2, this implies that the 
method is not SSP. □ 

We remark here that other results on the structure of SSP TSRK methods may 
be obtained by similar use of the stage order conditions and Theorems 2 and 3. We 
list some examples here, but omit the proofs since these results are not essential to 
our present purpose. 
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1. Any SSP TSRK method (implicit or explicit) of order greater than four must 
have a stage equal to it" -1 or u n . 

2. The abscissae c, of any SSP TSRK method of order greater than four must 
each be non-negative or equal to —1. 

3. (Implicit) SSP TSRK methods with p > 8 must be of Type II. 

4. Optimal SSP Two-step Runge Kutta methods. Our objective in this 
section is to find SSP TSRK methods that have the largest possible SSP coefficient. 
A method of order p with s stages is said to be optimal if it has the largest value of 
C over all TSRK methods with order at least p with no more than s stages. 

The methods presented were found via numerical search using Matlab's Opti- 
mization and Global Optimization toolboxes. We searched over Type I and Type II 
methods; the optimal methods found are of Type II in every case. For the methods of 
seventh and eighth order, this is known a priori from Theorem 4. Even for the lower 
order methods, this is not surprising, since explicit s-stage Type II methods (2.3) have 
an additional s — 1 degrees of freedom compared to explicit s-stage Type I methods 
(2.2). Although we do not know in general if these methods are globally optimal, our 
search recovered the global optimum in every case for which it was already known. 

4.1. Formulating the Optimization Problem. The optimization problem is 
formulated using the theory of Section 2: 

maxr, 

S,T 

r (I + rT)-^ > 0, 
subject to < (I + rT) _1 T > 0, 
{ $ p (S,T)=0, 

where the inequalities are understood component-wise and $p(S,T) represents the 
order conditions up to order p. This formulation, solved numerically in Matlab using 
a sequential quadratic programming approach (fmincon in the optimization toolbox), 
was used to find the methods given below. 

In comparing methods with different numbers of stages, one is usually interested 
in the relative time advancement per computational cost. For this purpose, we define 
the effective SSP coefficient 

s 

This normalization enables us to compare the cost of integration up to a given time, 
assuming that the time step is chosen according to the SSP restriction. 

Remark 4. By virtue of Theorem 1, optimal SSP methods found in the classes of 
Type I and Type II TSRK methods are in fact optimal over the larger class of methods 
(2.1). Also, because they do not use intermediate stages from previous timesteps, 
special conditions on the starting method (important for methods of the form (2.1) 
[16, 45, 44]) are unnecessary. Instead, the method can be started with any SSP Runge- 
Kutta method of the appropriate order. 

Remark 5. The optimal Type II methods found here could be rewritten as Type I 
methods by adding a stage (see Remark 1). However, this would reduce their effective 
SSP coefficient and render them (apparently) less efficient than some other (Type I) 
methods. This indicates once more the importance of explicitly accounting for the 
practical difference between Type I and Type II methods. 
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4.2. Low-storage implementation of Type II SSP TSRKs. The form 
(2.8), with r = C(S,T), typically yields very sparse coefficient matrices for optimal 
Type II SSP TSRK methods. This form is useful for a low-storage implementation. 
Written out explicitly, this form is: 



y ? = dV 1 - 1 + [ i - di - E m I «" + E (W + ' (1 - J - s) ' 

j=o / j=o V r / 

(4.1a) 

9u n -' , [ 1 - - E ^ ] «" + E ^ f»i + V^^^) ' (4Ab) 



j=0 / j=0 



where the coefficients are given by (using the relations (2.9)): 

Q = rA(I + rA)-\ n = rb T (I + rA)~\ 

d = d Qd, = -r| T d. 



When implemented in this form, many of the methods presented in the next 
section have modest storage requirements, despite using large numbers of stages. The 
analogous form for Runge-Kutta methods was used in [28]. 

In the following sections we discuss the numerically optimal methods, and in 
Tables A. 2, A. 3, A. 4, and A. 5, we give the coefficients in the low-storage form (4.1) 
for some numerically optimal methods. 

4.3. Optimal Methods of Orders One to Four. In the case of first-order 
methods, one can do no better (in terms of effective SSP coefficient) than the forward 
Euler method. For orders two to four, SSP coefficients of optimal methods found by 
numerical search are listed in Table 4.1. We list these mainly for completeness, since 
SSP Runge-Kutta methods with good properties exist up to order four. 

In [29], upper bounds for the values in Table 4.1 are found by computing optimally 
contractive general linear methods for linear systems of ODEs. Comparing the present 
results to the two-step results from that work, we see that this upper bound is achieved 
(as expected) for all first and second order methods, and even for the two- and three- 
stage third-order methods. 

Optimal methods found in [8] include two-step general linear methods of up to 
fourth order using up to four stages. By comparing Table 4.1 with the results therein, 
we see that the SSP coefficients of the optimal methods among the classes examined 
in both works (namely, for 1 < s < 4, 2 < p < 4) agree. The methods found in [8] are 
produced by software that guarantees global optimality. 

All results listed in bold are thus known to be optimal because they match those 
obtained in [8], [31], or both. This demonstrates that our numerical optimization 
approach was able to recover all known globally optimal methods, and suggests that 
the remaining methods found in the present work may also be globally optimal. 

The optimal s-stage, second-order SSP TSRK method is in fact both a Type I and 
Type II method, and was found in numerical searches over methods of both forms. It 
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Table 4.1 

Effective SSP coefficients C e ff of optimal explicit 2-step Runge-Kutta methods of order two to 
four. Results known to be optimal from [8] or [31] are shown in bold. 



s\p 


2 


3 


4 


2 


0.707 


0.366 




3 


0.816 


0.550 


0.286 


4 


0.866 


0.578 


0.398 


5 


0.894 


0.598 


0.472 


6 


0.913 


0.630 


0.509 


7 


0.926 


0.641 


0.534 


8 


0.935 


0.653 


0.562 


9 


0.943 


0.667 


0.586 


10 


0.949 


0.683 


0.610 



has SSP coefficient C = -J s(s — 1) and nonzero coefficients 

q%,i-i = 1, (2 < i < s), 

q s +i,s = 2(C - s+ 1), 
d = 0, 

9 = 2(s-C) - 1. 

Note that these methods have C e ff = y^^j whereas the corresponding optimal 

Runge-Kutta methods have C e ff = s= ^- Using the low-storage assumption intro- 
duced in [28], these methods can be implemented with just three storage registers, 
just one register more than is required for the optimal second-order SSP Runge-Kutta 
methods. 

The optimal nine-stage, third-order method is remarkable in that it is a Runge- 
Kutta method. In other words, allowing the freedom of using an additional step does 
not improve the SSP coefficient in this case. 

4.4. Optimal Methods of Orders Five to Eight. Table 4.2 lists effective SSP 
coefficients of numerically optimal TSRK methods of orders five to eight. Although 
these methods require many stages, it should be remembered that high-order (non- 
SSP) Runge-Kutta methods also require many stages. Indeed, some of our SSP 
TSRK methods have fewer stages than the minimum number required to achieve the 
corresponding order for an Runge-Kutta method (regardless of SSP considerations). 

The fifth-order methods present an unusual phenomenon: when the number of 
stages is allowed to be greater than eight, it is not possible to achieve a larger effective 
SSP coefficient than the optimal 8-stage method, even allowing as many as twelve 
stages. This appears to be accurate, and not simply due to failure of the numerical 
optimizer, since in the nine-stage case the optimization scheme recovers the apparently 
optimal method in less than one minute, but fails to find a better result after several 
hours. 

The only existing SSP methods of order greater than four are the hybrid methods 
of Huang [21]. Comparing the best TSRK methods of each order with the best 
hybrid methods of each order, the TSRK methods have substantially larger effective 
SSP coefficients. 

The effective SSP coefficient is a fair metric for comparison between methods 
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Table 4.2 

Effective SSP coefficients C e ff of optimal explicit two-step Runge-Kutta methods of order five 
to eight. 



s\p 


5 


6 


7 


8 


4 


0.214 








5 


0.324 








6 


0.385 


0.099 






7 


0.418 


0.182 






8 


0.447 


0.242 


0.071 




9 


0.438 


0.287 


0.124 




10 


0.425 


0.320 


0.179 




11 


0.431 


0.338 


0.218 


0.031 


12 


0.439 


0.365 


0.231 


0.078 



of the same order of accuracy. Furthermore, our twelve-stage TSRK methods have 
sparse coefficient matrices and can be implemented in the low-storage form (4.1). 
Specifically, the fifth- through eighth-order methods of twelve stages require only 
5, 7, 7, and 10 memory locations per unknown, respectively, under the low-storage 
assumption employed in [28, 30]. Typically the methods with fewer stages require the 
same or more storage, so there is no reason to prefer methods with fewer stages if 
they have lower effective SSP coefficients. Thus, for sixth through eighth order, the 
twelve-stage methods seem preferable. The SSP TSRK methods recommended here 
even require less storage than what (non-SSP one-step) Runge-Kutta methods of the 
corresponding order would typically use. 

In the case of fifth order methods, the eight-stage method has a larger effective 
SSP coefficient than the twelve stage method, so the eight stage method seems best in 
terms of efficiency. However the eight stage method requires more storage registers (6) 
than the twelve stage method (5) . So while the eight stage method might be preferred 
for efficiency, the twelve stage method is preferred for low storage considerations. 

5. Numerical Experiments. 

5.1. Start-up procedure. As mentioned in Section 4.1, TSRK are not self- 
starting and thus require startup procedures, and while in general this is somewhat 
complicated, for our class of methods it is straightforward. We only require that 
the starting procedure be of sufficient accuracy and that it also be strong stability 
preserving. 

Figure 5.1 demonstrates one possible start-up procedure that we employed in our 
convergence studies and our other numerical tests to follow. The first step of size At 
from to to t\ is subdivided into substcps in powers of two. The SSPRK(5,4) scheme 
[42, 32] or the SSPRK(10,4) scheme [28] is used for the first substep, with the stepsize 
At* chosen small enough so that the local truncation error of the Runge-Kutta scheme 
is smaller than the global error of the TSRK scheme. Specifically, this can be achieved 
for an TSRK method of order p = 5, 6, 7 or 8 by taking 

Af = ^, 7 6 Z, and (At*) 5 = AAt p = 0(At p ). (5.1) 

Subsequent substeps are taken with the TSRK scheme itself, doubling the stepsizes 
until reaching t\. From there, the TSRK scheme repeatedly advances the solution 
from t n to t n+ \ using previous step values u n -\ and u n . 
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Fig. 5.1. One possible startup procedure for SSP TSRK schemes. The first step from to to t± is 
subdivided into substeps (here there are three substeps of sizes ^, and ^). An SSP Runge-Kutta 
scheme is used for the first substep. Subsequent substeps are taken with the TSRK scheme itself, 
doubling the stepsizes until reaching t\ . We emphasize that the startup procedure is not critical for 
this class of TSRK methods. 




Fig. 5.2. Convergence results for some TSRK schemes on the Dahlquist test problem (left) and 
van der Pol problem (right). The slopes of the lines confirm the design orders of the TSRK methods. 



5.2. Order Verification. Convergence studies on two ODE test problems con- 
firm that the SSP TSRK methods achieve their design orders. The first is the 
Dahlquist test problem u' = \u, with u = 1 and A = 2, solved until tt = 1. 
Figure 5.2 shows a sample of TSRK methods achieving their design orders on this 
problem. The starting procedure used SSPRK(fO,4) with the constant A in (5.f ) set 
respectively to [~, ~, 1CT 2 , 1CT 3 , 10~ 3 ] for orders p = 4, 5, 6, 7, and 8. 

The nonlinear van der Pol problem (e.g., [35]) can be written as an ODE initial 
value problem consisting of two components 

u[ = u 2 , 

U '-2 = \ + ( 1 _ U l) M 2) , 

where we use e = 0.01 with corresponding initial conditions u° = [2; —0.6654321] and 
solve until tf — \- The starting procedure used SSPRK(10,4) with constant A = 1 in 
(5.1). Error in the maximum norm is estimated against a highly-accurate reference 
solution calculated with Matlab's ODE45 routine. Figure 5.2 shows a sample of the 
TSRK schemes achieving their design orders on this problem. 

5.3. High-order WENO. Weighted essentially non-oscillatory schemes (WENO) 
[18, 17, 26] are finite difference or finite volume schemes that use linear combination 
of lower order fluxes to obtain a higher order approximations, while ensuring non- 
oscillatory solutions. This is accomplished by using adaptive stencils which approach 
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Fig. 5.3. Convergence results for 2D advection using rth-order WENO discretizations and the 
TSRK integrators (c.f., [13, Fig. 15]). Maximum error versus number of spatial grid points in each 
direction (left). Observed orders of accuracy calculated from these errors (right). Computed using 
the same parameters as [13, Fig. 15] (final time tf = 20, At = 0.5 Ax, mapped WENO spatial 
discretization with pp = r). Starting procedure as described in Section 5.1 using the SSPRK(5,4) 
scheme for the initial substep. 



centered difference stencils in smooth regions and one-sided stencils near discontinu- 
ities. Many WENO methods exist, and the difference between them is in the compu- 
tation of the stencil weights. WENO methods can be constructed to be high order 
[13, 3]. In [13], WENO of up to 17th-order were constructed and tested numerically. 
However, the authors note that in some of their computations the error was limited by 
the order of the time integration, which was relatively low (third-order SSPRK(3,3)). 
In Figure 5.3, we reproduce the numerical experiment of [13, Fig. 15], specifically the 
2D linear advection of a sinusoidal initial condition uq(x, y) — sin(7r(x + y)), in a 
periodic square using various high-order WENO methods and our TSRK integrators 
of order 5, 7 and 8 using 12 stages. Compared with [13, Fig. 15], we note that the 
error is no longer dominated by the temporal error. Thus the higher-order SSP TSRK 
schemes allow us to see the behavior of the high-order WENO spatial discretization 
schemes. 

5.4. Buckley Lever ett. The Buckley-Leverett equation is a model for two- 
phase flow through porous media and consists of the conservation law 

U t +f(U) x =0, with f(U)= u2 + ^_ ur 

We use a = h and initial conditions 

u(x 0) = i 1 iix ^h 
[0 otherwise, 

on x <E [0, 1) with periodic boundary conditions. Our spatial discretization uses 100 
points and following [23, 31] we use a conservative scheme with Koren limiter. We 
compute the solution until tf = g. For this problem, the Euler solution is total 
variation diminishing (TVD) for At < A£fe = 0.0025 [31]. As discussed above, we 
must also satisfy the SSP time-step restriction for the starting method. 

Figure 5.4 shows typical solutions using an TSRK scheme with timestep At = 
erAipE- Table 5.1 shows the maximal TVD time-step sizes, expressed as At = 
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t = 0.16625, TV = 1 t = 0.168, TV = 1.0306 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 



x x 

Fig. 5.4. Two numerical solutions of the Buckley-Leverett test problem. Left: time-step satisfies 
the SSP time-step restriction (TSRK(8,5) using At = 3.5AtpEj. Right: time-step does not satisfy 
the restriction (At = 5.6AtpEj and visible oscillations have formed, increasing the total variation 
of the solution. 

Table 5.1 

SSP coefficients versus largest time steps exhibiting the TVD property (At = cr BL Atp^) on the 
Buckley-Leverett example, for some of the SSP TSRK(s,p) schemes. The effective SSP coefficient 
C e ff should be a lower bound for crg^/s and indeed this is observed. SSPRK(10,4) [28] is used as 
the first step in the starting procedure. 



Method 


theoretical 


observed 




C 


C c s 


0"BL 


o"bl/s 


TSRK(4,4) 


1.5917 


0.398 


2.16 


0.540 


TSRK(8,5) 


3.5794 


0.447 


4.41 


0.551 


TSRK(12,5) 


5.2675 


0.439 


6.97 


0.581 


TSRK(12,6) 


4.3838 


0.365 


6.80 


0.567 


TSRK(12,7) 


2.7659 


0.231 


4.86 


0.405 


TSRK(12,8) 


0.94155 


0.0785 


4.42 


0.368 



ctbl A£fe, for the Buckley-Leverett test problem. The results show that the SSP 
coefficient is a lower bound for what is observed in practice, confirming the theoreti- 
cal importance of the SSP coefficient. 

6. Conclusions. In this paper we have analyzed the strong stability preserv- 
ing property of two-step Runge-Kutta (TSRK) methods. We find that SSP TSRK 
methods have a relatively simple form and that explicit methods are subject to a 
maximal order of eight. We have presented numerically optimal SSP TSRK meth- 
ods of order up to this bound of eight. These methods overcome the fourth order 
barrier for (one-step) SSP Runge-Kutta methods and allow larger SSP coefficients 
than the corresponding order multi-step methods. The discovery of these methods 
was facilitated by our formulation of the optimization problem in an efficient form, 
aided by simplified order conditions and constraints on the coefficients derived by 
using the SSP theory for general linear methods. These methods feature favorable 
storage properties and are easy to implement and start up, as they do not use stage 
values from previous steps. 

We show that high-order SSP two-step Runge-Kutta methods are useful for the 
time integration of a variety of hyperbolic PDEs, especially in conjunction with high- 
order spatial discretizations. In the case of a Buckley-Leverett numerical test case, 
the SSP coefficient of these methods is confirmed to provide a lower bound for the 
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actual time-step needed to preserve the total variation diminishing property. 

The order conditions and SSP conditions we have derived for these methods ex- 
tend in a very simple way to methods with more steps. Future work will investigate 
methods with more steps and will further investigate the use of start-up methods for 
use with SSP multi-step Runge-Kutta methods. 

Acknowledgment. The authors are grateful to an anonymous referee, whose careful 
reading and detailed comments improved several technical details of the paper. 

Appendix A. Coefficients of Numerically Optimal Methods. 

Table A.l 

Coefficients of the optimal explicit 8-stage 5th-order SSP TSRK method (Type II) 



8 = 

d = 1.000000000000000 
d 7 = 0.003674184820260 
r)2 = 0.179502832154858 
7j3 = 0.073789956884809 
7j 6 = 0.017607159013167 
rjs = 0.729100051947166 
q 2 o = 0.085330772947643 
q 3 = 0.058121281984411 



97 a = 0.020705281786630 
qs[o = 0.008506650138784 
q 2il = 0.914669227052357 
94,1 = 0.036365639242841 
q bl = 0.491214340660555 
q 6A = 0.566135231631241 
9 7|1 = 0.091646079651566 

9 8 i = 0.110261531523242 
q 3 2 = 0.941878718015589 



q s 2 = 0.030113037742445 
q 4 ,3 = 0.802870131352638 
q 5A = 0.508785659339445 
9 6j5 = 0.433864768368758 
97j6 = 0.883974453741544 
98,7 = 0.851118780595529 



Table A. 2 

Coefficients of the optimal explicit 12-stage 5th-order SSP TSRK method (Type II) 



= 
do = 1 

rji = 0.010869478269914 
Tie = 0.252584630617780 
7710 = 0.328029300816831 
T/i2 = 0.408516590295475 
g 2j0 = 0.037442206073461 
q 3 \o = 0.004990369159650 
q 2A = 0.962557793926539 



q 6 1 = 0.041456384663457 
q 7 \i = 0.893102584263455 
99,1 = 0.103110842229401 
q 10A = 0.109219062395598 
qil ' A = 0.069771767766966 
q 12 ' tl = 0.050213434903531 
q 3t 2 = 0.750941165462252 
q 4 3 = 0.816192058725826 
95,4 = 0.881400968167496 



q 6i5 = 0.897622496599848 
q 7 [ e = 0.106897415736545 
98,6 = 0.197331844351083 
q SJ = 0.748110262498258 
q 9 ] 8 = 0.864072067200705 
q 10 9 = 0.890780937604403 
9ii]io = 0.928630488244921 
912,11 = 0.949786565096469 



Table A. 3 

Coefficients of the optimal explicit 12-stage 6th-order SSP TSRK method (Type II) 



8 = 2. 455884612148108c - 04 
do = 1 

diO = 0.000534877909816 
92 o = 0.030262100443273 
92^1 = 0.664746114331100 
96i = 0.656374628865518 
97 1 = 0.210836921275170 
99^! = 0.066235890301163 
910,1 = 0.076611491217295 
g 12j i = 0.016496364995214 
9 32 = 0.590319496200531 



94,3 = 0.729376762034313 
95^4 = 0.826687833242084 

910.4 = 0.091956261008213 
9n j4 = 0.135742974049075 
96 5 = 0.267480130553594 

911.5 = 0.269086406273540 
9 12 5 = 0.344231433411227 
9 7>6 = 0.650991182223416 

912.6 = 0.017516154376138 
98,7 = 0.873267220579217 
99^8 = 0.877348047199139 



910.9 = 0.822483564557728 

911.10 = 0.587217894186976 

912.11 = 0.621756047217421 
Tji = 0.012523410805564 

rjs = 0.094203091821030 
jj 9 = 0.318700620499891 
rjio = 0.107955864652328 
tj 12 = 0.456039783326905 
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Table A. 4 

Coefficients of the optimal explicit 12-stage 7th-order SSP 



TSRK method (Type II) 



8 = 1. 040248277612947c - 04 
d = 1.000000000000000 
d 2 = 0.003229110378701 
d 4 = 0.006337974349692 
d 5 = 0.002497954201566 
lis = 0.017328228771149 
di2 = 0.000520256250682 
m = 0.000515717568412 
tji = 0.040472655980253 
Ve = 0.081167924336040 
r; 7 = 0.238308176460039 
r te = 0.032690786323542 
n 12 = 0.547467490509490 
q 2 o = 0.147321824258074 



q 2 1 = 0.849449065363225 
q 3 [i = 0.120943274105256 
g 4 ,i = 0.368587879161520 
g 5 ,i = 0.222052624372191 
q 6A = 0.137403913798966 
= 0.146278214690851 
98 ! = 0.444640119039330 
q 9 1 = 0.143808624107155 
(jio.i = 0.102844296820036 
gn,i = 0.071911085489036 
912,1 = 0.057306282668522 
93, 2 = 0.433019948758255 
q 7 \ 2 = 0.014863996841828 
99,2 = 0.026942009774408 



94.3 = 0.166320497215237 
9i ,3 = 0.032851385162085 

95.4 = 0.34 3 703 78075 9466 
q 6 [ 5 = 0.519758489994316 
q 7 [ e = 0.598177722195673 
q s 'j = 0.488244475584515 
910 7 = 0.356898323452469 
9n, 7 = 0.508453150788232 
912,7 = 0.496859299069734 
99, 8 = 0.704865150213419 
910,9 = 0.409241038172241 
?n,io = 0.327005955932695 
912,11 = 0.364647377606582 



Table A. 5 

Coefficients of the optimal explicit 12-stage 8th-order SSP TSRK method (Type II) 



8 = 4. 796147528566197c - 05 
do = 1.000000000000000 
d 2 = 0.036513886685777 
d 4 = 0.004205435886220 
d & = 0.000457751617285 
d 7 = 0.007407526543898 
d 8 = 0.000486094553850 
Tji = 0.033190060418244 
rj2 = 0.001567085177702 
?J3 = 0.014033053074861 
Vi = 0.017979737866822 
775 = 0.094582502432986 
r/6 = 0.082918042281378 
r; 7 = 0.020622633348484 
7 )8 = 0.033521998905243 
7(9 = 0.092066893962539 
T/io = 0.076089630105122 
7)ii = 0.070505470986376 
77i2 = 0.072975312278165 



92 = 0.017683145596548 
9 3 ^ = 0.001154189099465 
9 6 ^ = 0.000065395819685 

99.0 = 0.000042696255773 
9^,0 = 0.000116117869841 
912^0 = 0.000019430720566 

92.1 = 0.154785324942633 
94,1 = 0.113729301017461 
9 5 'i = 0.061188134340758 
96^1 = 0.068824803789446 

97.1 = 0.133098034326412 
9 8 ^i = 0.080582670156691 
99^! = 0.038242841051944 
910,1 = 0.071728403470890 
911,1 = 0.053869626312442 
9i 2 ,i = 0.009079504342639 

93.2 = 0.200161251441789 
9 6 ' 2 = 0.008642531617482 
9 4 ' 3 = 0.057780552515458 



99.3 = 0.029907847389714 

95.4 = 0.165254103192244 
97^4 = 0.005039627904425 
98,4 = 0.069726774932478 

99.4 = 0.022904196667572 
912,4 = 0.130730221736770 

96.5 = 0.229847794524568 
99^ 5 = 0.095367316002296 
976 = 0.252990567222936 
99^6 = 0.176462398918299 
910 6 = 0.281349762794588 
911,6 = 0.327578464731509 
912,6 = 0.149446805276484 

98.7 = 0.324486261336648 

99.8 = 0.120659479468128 
9io,g = 0.166819833904944 
9n,i = 0.157699899495506 
912,11 = 0.314802533082027 
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